to remove/masterdosresNMA.R

# libraries
library(devtools)
install_github("htx-r/doseresNMA",force=TRUE)
library(doseresNMA)

library(R2jags)
library(rms)
library(meta)
library('MBNMAdose')

#** Prepare data

# load
mydata <- read.csv('DOSEmainanalysis.csv')
antidepdata0=mydata[mydata$exc==F,]
antidepdata <- antidepdata0[antidepdata0$Study_No!=161001,]
length(table(antidepdata0$Dose_delivered_mean))
# add OR to the dataset
antidepdata$studyid <- as.numeric(as.factor(antidepdata$Study_No))
antidepdata$nonResponders <- antidepdata$No_randomised- antidepdata$Responders
logORmat <- sapply(unique(antidepdata$studyid),function(i) createORreference.fun(antidepdata$Responders[antidepdata$studyid==i],antidepdata$No_randomised[antidepdata$studyid==i]),simplify = FALSE)
logORmat <- do.call(rbind,logORmat)
antidepdata$logOR <- c(logORmat[,1])
antidepdata$selogOR <- c(logORmat[,2])

#** JAGS data format
class <- as.numeric(as.factor(antidepdata$Drug))
class[class==5] <- 'A'
class[class==1|class==2|class==3] <- 'B'
class[class==4|class==6] <- 'C'
class <- as.numeric(as.factor(class))
mydata <- data.frame(studyid=antidepdata$Study_No,drug=antidepdata$Drug,dose=antidepdata$Dose_delivered_mean,
                   r=antidepdata$Responders,n=antidepdata$No_randomised)
                   # ROB=antidepdata$Overall_study_RoB,Study_year=antidepdata$Study_year,selogOR=antidepdata$selogOR,class=class)
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=F,class=F,ref='placebo')

# knots=c(20,30,50)
# dosresnetmeta.jagsdata$nn <- dosresnetmeta.jagsdata$n
# dosresnetmeta.jagsdata$rr <- dosresnetmeta.jagsdata$r
# dosresnetmeta.jagsdata$new.dose <-  seq(1,80,1)
# dosresnetmeta.jagsdata$f.new.dose <- rcspline.eval(dosresnetmeta.jagsdata$new.dose,knots,inclx = T)[,2]
# dosresnetmeta.jagsdata$nd.new <- length(dosresnetmeta.jagsdata$new.dose)

#** Run JAGS

# 1. Run DRNMA - cubic splines
#dosresnetmeta.jagsdata$ref <- 15
# dosresnetmeta.jagsdata$trt1 <- unique(dosresnetmeta.jagsdata$t[,1])
# dosresnetmeta.jagsdata$trt2 <- unique(c(dosresnetmeta.jagsdata$t[,-1]))[!is.na(unique(c(dosresnetmeta.jagsdata$t[,-1])))]
dosresnetmeta.jagsdata$nd.new
dosresnetmeta.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','dd','dd2','tau','Z','p.drug','cumeffectiveness'),model.file = modelDRnetmetaSpline,
                                          n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmeta.runjagsRCSInc <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('dd','dd2','tau','Z','p.drug','cumeffectiveness'),model.file = modelDRnetmetaSplineInc,
                                          n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmeta.runjagsRCS$BUGSoutput$summary
names.inc1 <- paste0('dd[',rep(dosresnetmeta.jagsdata$trt1,each=length(dosresnetmeta.jagsdata$trt2)),',',dosresnetmeta.jagsdata$trt2,']')
names.inc2 <- paste0('dd2[',rep(dosresnetmeta.jagsdata$trt1,each=length(dosresnetmeta.jagsdata$trt2)),',',dosresnetmeta.jagsdata$trt2,']')

est.con1 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[rownames(dosresnetmeta.runjagsRCS$BUGSoutput$summary)%in%names.inc1,]
est.con2 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[rownames(dosresnetmeta.runjagsRCS$BUGSoutput$summary)%in%names.inc2,]

est.incon1 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[names.inc1,]
est.incon2 <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[names.inc2,]

plot(abs(est.con1[,'mean']),abs(est.incon1[,'mean']),ylim = c(0,0.03),xlim = c(0,0.03))
abline(0,1)
myd <- data.frame(y=est.con2[,'mean'],x=est.incon2[,'mean'])
ggplot(myd,aes(y=y,x=x))+
  geom_point()+
  geom_abline(slope=1,intercept=0)+
  theme_bw()
# 2. Run DRNMR - cubic splines - ROB
mydata$covariate <- as.numeric(as.factor(antidepdata$Overall_study_RoB))-2
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=T)
dosresnetmeta.jagsdata$ref <- 5

dosresnetmetaregROB.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','bb','tua'),model.file = modelDRnetmetaregSpline,
                                             n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregROB.runjagsRCS

# 3. Run DRNMR - cubic splines - VAR: small-study effect
mydata$covariate <- antidepdata$selogOR^2
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=T)
dosresnetmeta.jagsdata$ref <- 5

dosresnetmetaregVAR.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','bb','tau'),model.file = modelDRnetmetaregSpline,
                                                n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregVAR.runjagsRCS

# 4. Run DRNMR - cubic splines - year of study **** Does not run, study_year has NA's
mydata$covariate <- antidepdata$Study_year-min(antidepdata$Study_year,na.rm=T)
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata[!is.na(mydata$covariate),],cov=T)
dosresnetmeta.jagsdata$ref <- 5

dosresnetmetaregYEAR.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','bb','tau'),model.file = modelDRnetmetaregSpline,
                                                n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregYEAR.runjagsRCS


# 5. Run DRNMA - cubic splines + class effect
mydata$class <- class
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=F,includeClass=T)
dosresnetmeta.jagsdata$ref <- 5

dosresnetmetaClass.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('dd','dd2','sd'),model.file = modelDRnetmetaSplineClass,
                                                 n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaClass.runjagsRCS

# 6. Run DRNMR - cubic splines + class effect + covariate - ROB
mydata$covariate <- as.numeric(as.factor(antidepdata$Overall_study_RoB))-2
dosresnetmeta.jagsdata <- makejagsDRnetmeta(data=mydata,cov=T,includeClass=T)
dosresnetmeta.jagsdata$nc <- 3
dosresnetmeta.jagsdata$ref <- 5

dosresnetmetaregClass.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('dd','dd2','sd'),model.file = modelDRnetmetaregSplineClass,
                                                  n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
dosresnetmetaregClass.runjagsRCS

#
as.numeric(as.factor(antidepdata$Drug))
#
sort(levels(as.factor(antidepdata$Drug)))
#
sort(c('c','b','a'))

install_github("audrey-b/BUGSnet",force=TRUE)
library(BUGSnet)
data("thrombolytic")
head(thrombolytic)
myd <- data.prep(antidep,varname.t = 'drug',varname.s = 'studyid')
net.plot(myd,node.scale = 0.2,edge.scale = 0.25)
#


# Fractional polynomials
library(mfp)
library(dosresmeta)
str(antidepdata)
f <- mfp(logRR~fp(hayasaka_ddd),data = antidepdata,family = gaussian)
summary(f)
plot(f)
antidepdata$hayasaka_ddd2 <- ifelse(antidepdata$hayasaka_ddd ==0, 0,antidepdata$hayasaka_ddd^-3)
ff <- dosresmeta(formula=logRR~I(hayasaka_ddd)+hayasaka_ddd2, proc="1stage",id=Study_No, type=type,cases=Responders,n=No_randomised,se=selogRR,data=antidepdata,method = 'reml')
ff2 <- dosresmeta(formula=logRR~rcs(hayasaka_ddd,knots), proc="1stage",id=Study_No, type=type,cases=Responders,n=No_randomised,se=selogRR,data=antidepdata,method = 'reml')

#
knots=c(20,30,50)
newdose <- 1:80
dd <- 0.0058*newdose+10.4571*newdose^-2
dd1 <- 65.8204*newdose^-2
dd3 <- 0.0001*newdose^2+46.9503*newdose^-2
dd2 <- 0.0091*newdose+-0.0156*rcspline.eval(newdose,knots,inclx = T)[,2]
with(antidepdata,plot(hayasaka_ddd,logRR))
with(antidepdata,lines(dd))
with(antidepdata,lines(dd1,col=3))
with(antidepdata,lines(dd2,col=2))
with(antidepdata,lines(dd3,col=4))

#








# predict the absolute effect + ranking: SUCRA

# dosresnetmeta.runjagsRCS <- jags.parallel(data = dosresnetmeta.jagsdata,inits=NULL,parameters.to.save = c('d','d2','Z','sd','p.drug','SUCRA'),model.file = modelDRnetmetaSpline,
#                                           n.chains=3,n.iter = 1000,n.burnin = 200,DIC=F,n.thin = 1)
# dosresnetmeta.runjagsRCS$BUGSoutput$mean$SUCRA
#
# p.placebo<- exp(dosresnetmeta.runjagsRCS$BUGSoutput$mean$Z)/(1+exp(dosresnetmeta.runjagsRCS$BUGSoutput$mean$Z))
# p.drug <- dosresnetmeta.runjagsRCS$BUGSoutput$mean$p.drug
# l.ci <-  dosresnetmeta.runjagsRCS$BUGSoutput$summary[paste0('p.drug[',1:80,',',6,']'),'2.5%']
# u.ci <- dosresnetmeta.runjagsRCS$BUGSoutput$summary[paste0('p.drug[',1:80,',',6,']'),'97.5%']
# lp.ci <-  exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','2.5%'])/(1+exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','2.5%']))
# up.ci <- exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','97.5%'])/(1+exp(dosresnetmeta.runjagsRCS$BUGSoutput$summary['Z','97.5%']))
#
# library(ggplot2)
# df1 <- data.frame(new.dose=dosresnetmeta.jagsdata$new.dose,y1 =p.drug[,5],y2=l.ci,
#                   y3=u.ci,yp1=p.placebo,yp2=lp.ci,yp3=up.ci)
#
# ggplot(data = df1,aes(x=new.dose)) +
#   geom_line(aes(y=y1,color='treatment'))+
#   geom_line(aes(y=yp1,color='placebo'))+
#   scale_color_manual(values=c('steelblue','coral4'))+
#   geom_smooth(aes(x=new.dose, y=y1, ymin=y2,
#                   ymax=y3),color='coral4',fill='coral1',
#               data=df1, stat="identity")+
#   geom_smooth(aes(x=new.dose, y=yp1, ymin=yp2,
#                   ymax=yp3),color='steelblue',fill='steelblue2',
#               data=df1, stat="identity")+
#   xlab('')+
#   ylab('')+
#   ylim(0,0.9)+
#   theme(panel.background = element_rect(fill = 'snow1',colour = 'white'),
#         legend.position = c(0.72,0.95), legend.key.size = unit(3, "cm"), legend.text = element_text(size=12),
#         legend.key.height  = unit(0.5,"cm"),legend.title = element_blank(), legend.background = element_blank(),
#         axis.text.x = element_text(face='bold',size=14),
#         axis.text.y = element_text(face='bold',size=14))









# Timesheet_Template_2019 <- read_excel("~/Google Drive/Administrative/Timesheet_Template_2019.xlsx",
#                                         sheet = "JAN", col_names = FALSE)
# Timesheet_Template_2019 <- as.data.frame(Timesheet_Template_2019)
# Timesheet_Template_2019[c('Work Package 1'),2:(ncol(Timesheet_Template_2019)-2)] <- sample(seq(7,9,by =0.5),ncol(Timesheet_Template_2019)-3,replace = T)


# for (i in dosresnetmeta.jagsdata$trt2) {
#   print(i)
# }
#
htx-r/doseresNMA documentation built on Jan. 28, 2021, 5:32 a.m.